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Rayleigh functional for nonlinear systems 

Valery S. ShchesnoviciQ and Solange B. Cavalcanti 
Departamento de Fisica - Universidade Federal de Alagoas, Maceid AL 57072-970, Brazil 

Abstract 

We introduce Rayleigh functional for nonlinear systems. It is defined using the energy func- 
tional and the normalization properties of the variables of variation. The key property of the 
Rayleigh quotient for linear systems is preserved in our definition: the extremals of the Rayleigh 
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functional coincide with the stationary solutions of the Euler-Lagrange equation. Moreover, the 
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second variation of the Rayleigh functional defines stability of the solution. This gives rise to a 
powerful numerical optimization method in the search for the energy minimizers. It is shown that 
the well-known imaginary time relaxation is a special case of our method. To illustrate the method 
we find the stationary states of Bose-Einstein condensates in various geometries. Finally, we show 
that the Rayleigh functional also provides a simple way to derive analytical identities satisfied by 
the stationary solutions of the critical nonlinear equations. 
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I. INTRODUCTION 



In dealing with systems of equations having a large number of degrees of freedom one 
frequently adopts an approximation which leads to a nonlinear partial differential equation. 
For instance, the mean-field approach in the statistical physics is based on the introduction of 
an order parameter governed by a nonlinear equation. For example, the mean-field theory of 
the Bose-Einstein condensate of a degenerate quantum gas is based on the Gross-Pit aevskii 
equation The Gross-Pitaevskii equation is the nonlinear Schrodinger equation with 
an external potential. It describes the so-called "matter waves" - the nonlinear collective 
modes of the degenerate quantum gas below the condensation temperature. In many other 
cases the approximate nonlinear equations, appearing as the leading order description of 
nonlinear waves in various branches of modern physics, for instance, in nonlinear optics, 
plasma, and on water, are of the nonlinear Schrodinger class (see, for instance, Ref. 2j and 
the references therein). In physics one is especially interested in the stationary solutions 
(a.k.a. the stationary points) of the governing equations, in particular, in the stationary 
points which minimize the energy. Only in some exceptional cases the solution can be 
obtained analytically and one has to rely on numerical simulations. If the nonlinear equation 
in question possesses the Lagrangian formulation then the stationary points can be found 
by a nonlinear optimization method. 

We propose a new optimization method for nonlinear partial differential equations, in 
particular, for the equations of the nonlinear Schrodinger (NLS) type. The method consists 
of a numerical minimization of the Rayleigh functional. In a broader perspective, we in- 
troduce the concept of Rayleigh functional for nonlinear systems and show how it can be 
employed in both numerical and analytical analysis of equations of physical importance. 

In our case, the optimization problem is to find the ground state of a nonlinear system, 
i.e. the stationary solution which minimizes the energy. We work in the space of smooth 
functions ip(x, t) of n spatial variables x = (x\, x n ) and time variable, which have bounded 
/2-norm: ||-?/>|| 2 = J d n x\ip\ 2 < 00. Thus we consider only bounded stationary solutions which 
decay to zero as Xk — > 00 and correspond to finite energies. 

We adopt the following energy functional: 

E{W} = J d n x£(x,ij(x,t),r(x,t),VtP(x,t),Vr(x,t)) (1) 
and assume that it has the phase invariance symmetry if) — ► e %e ip (we use the complex 



conjugate functions ip and tp* as independent variational variables instead of the real and 
imaginary parts of ip). The phase invariance symmetry results in conservation of the ^-norm. 

Our method also works for the generalization of the energy functional (£[]) to several 
variables of variation: tpk(x), k = 1, . . . ,m (this case is discussed below), and to the higher 
order derivatives: V p ip, p = 1, . . . ,s. We adopt the notations used in statistical physics, an 
important area for applications. Hence, the dependent variable ip(x ,t) will be referred to 
as the order parameter. 

The well known Gross-Pitaevskii functional for the order parameter of a Bose-Einstein 
condensate, 

E = J d 3 x ||^|V^(x,t)| 2 + V(f)|^(x,t)| 2 + ^|^(x,t)| 4 |, (2) 

belongs to the class specified by equation (JTJ). Here x = (x, y,z), N is the number of atoms 
in the condensate, g = 4irh 2 a s /m is the interaction coefficient due to the s-wave scattering 
of the atoms and V(x) is the trap potential (created by a magnetic field or non-resonant 
laser beams). The order parameter is normalized to 1 in equation (J2J). 

Let us briefly recall the basics of the nonlinear optimization. One is interested in the 
stationary state, given as ip(x,t) = e~ tflt ty(x), where \i is the chemical potential. The 
Euler-Lagrange equation corresponding to the energy functional (JTJ) reads 

,cty> _SE _d£ d£ 
1 dt 5ip* d^)* dVip* ' 

Here (and below) we denote 5F/5ijj* the variational derivative of a functional F with respect 
to ip* (thus we part with the usual notation of the latter by using the symbol V, while V 
is reserved for the usual gradient of a function). The stationary state satisfies the equation 
= SE/S^*. The idea of the imaginary time relaxation method is based on the fact 
that the variational derivative of a functional is an analog of the gradient of a function. 
Thus, by introducing the "imaginary time" r = it in equation (jHJ) one forces the order 
parameter to evolve in the direction of maximum decrease in the energy. The attractor 
of such evolution is, hopefully, a local minimum (in general, just a stationary point). The 
imaginary time evolution does not conserve the Z 2 - n orm and one must normalize the order 
parameter directly during the evolution. In other words, one allows the order parameter to 
evolve in the space of arbitrary functions but normalizes the solution after each step by the 
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prescription tp — > tAt. This can be formulated in a single equation: 

^ _ 5E{f,n 



dr Sf* 



(4) 



f_ V 



The method of Lagrange multipliers can be used to take into account the constraints 
imposed on the variables of variation (for instance, the fixed norm). It consists of a numerical 
minimization of an augmented functional which is the energy functional plus the constraint 
with a Lagrange multiplier. The combined functional has the stationary solution as its 
extremal point and the constrained minimization problem is converted into an unconstrained 
one. In this case, the time evolution of equation (jlj) can be substituted by a finite-step 
minimization scheme, such as the method of steepest descent or the conjugate gradient 
method, supplemented by an appropriate line-search algorithm (consult, for instance, Refs. 
[3|, |4j). For instance, the following two functionals are used 

F x ty,P} = V>*}-A* / d^l^l 2 , F 2 {i,,r} = Ety,P}+± (a - J d n x , (5) 

where the variable of variation is arbitrary (has arbitrary e2-norm). Indeed, the first of these 
functionals evidently has the stationary solutions as its extremals, while the second has the 
variation SE/Sip* — (H^H 2 — A)^. Setting \i — (A — ||^|| 2 ) we get the stationary point by 
equating the variation of F 2 to zero. The use of the functional F 2 in the search for the energy 
minimizers was advocated recently in Ref. [fj. The point is that the trivial solution ip — 0, 
being an extremal, frequently makes reaching other solutions difficult with the use of the 
functional F%, while functional F 2 is free from such a flaw. 

The simplest numerical realization of the minimization method is given by the steepest 
descent scheme: 

^ +1 = fc _ A ^%M, (6) 

where the parameter f3k is selected by an appropriate line-search algorithm. 

This paper is organized as follows. In section |H] we introduce the Rayleigh functional 
for the energy given by equation (P) and discuss its properties as a variational functional. 
The generalization to several order parameters is also considered. In section ITTT1 we present 
examples of the nonlinear optimization based on the Rayleigh functional. As an application 
of the Rayleigh functional in the analytical approach, in section IIVI we relate an identity 
satisfied by the stationary solutions of the critical NLS equations to the scale invariance 
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symmetry of the Rayleigh functional. Finally, in section the advantage of the numerical 
optimization method based on the Rayleigh functional is discussed. 



II. DEFINITION AND PROPERTIES OF THE RAYLEIGH FUNCTIONAL 

Before formulating our nonlinear optimization method, it would be instructive to recall 
how the eigenfunctions of a linear operator can be obtained numerically. Consider, for 
instance, the textbook problem of finding the eigenvalues and eigenfunctions of the Hamil- 
tonian operator H = — V 2 + V(x ) for a quantum particle in a potential well (we use H — 1 
and m = 1/2). This problem can be^ reformulated as an optimization problem by employing 
the well-known Rayleigh quotient 



6] 



= ^f^ff* > = /ftfif (7) 



J d 3 x \tp{x 

Here the function of variation ip is arbitrary, i.e. not normalized. The eigenfunctions are 
the stationary points of the Rayleigh quotient. Indeed 

|£ = pji (**-*) = "■ < 8 > 

where the eigenvalue is given as e = 7Z{ip, ip*}, i.e. it is the value of 1Z at the eigenfunction. 

Note that the Rayleigh quotient allows one to cast the constrained minimization problem 
(with the constraint ||"0|| = 1) into an unconstrained one. 



A. Rayleigh functional for a single order parameter 

Stationary solutions of nonlinear systems, the so-called nonlinear modes, can be consid- 
ered as the generalization of the eigenfunctions. They can be found in a way similar to the 
solution of the above eigenvalue problem. 

In contrast to the linear systems, a nonlinear one usually possesses continuous families 
of stationary solutions, each solution corresponding to a different chemical potential (an 
analog of the energy level in the nonlinear case). In general, the chemical potential takes 
its values from an infinite interval of the real line. By fixing the Z 2 - n orm we introduce it 
as a parameter in the energy functional. Thus, the chemical potential, being a continuous 
function of the norm, is also fixed (in general, we get away with a finite number of distinct 
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values lying on the different branches of the function /x = //(||^>||)). The only exception is 
the so-called critical case when a continuous family of stationary solutions corresponds to 
the same / 2 -norm (see also section HVI where we discuss the critical case). Therefore, apart 
from the critical case, by fixing the norm we select a finite number of solutions. Among 
them there is the energy minimizer (for a given value of the conserved / 2 -norm). We can 
reconstruct the continuous family of the stationary solutions by using different values of the 
norm. In the case of Bose-Einstein condensates, for instance, this corresponds to fixing the 
number of atoms and looking for the corresponding ground state solution. 

It is assumed that the nonlinear equation in question has a phase invariance resulting 
in the / 2 -norm conservation. In the following we will distinguish between the normalized 
and not-normalized functions, denoting the former by f(x,t) and the latter by ip(x,t). For 
convenience, we set the /2-norm of the stationary solution equal to 1, ||/|| = 1. This 
normalization is performed by using a scale transformation of the function of variation and 
results in the explicit appearance of the value of /2-norm in the energy functional (note that 
the chemical potential may also be scaled appropriately, as in equations (J2Hj) and (J2%)) of 
section ITl7|) . For example, in equation (j2J) we have the coefficient N at the nonlinear term, 
after we set the / 2 -norm of the solution to 1, and the value of the energy functional is, in 
fact, the energy of the condensate per atom. 

For the energy given by equation the following functional can serve as the nonlinear 
generalization of the Rayleigh quotient: 



Here the /2-norm | | is not fixed but is a functional of ip and ip* (since the norm is constant 
function of x we can take it out of the gradient operator V). Note that 1Z is a compound 
functional: it depends on the complex conjugate functions ip(x) and ip*(x) through the 
normalized ones, f{x) and f*(x), used in the energy functional E. This simple functional 
turns out to be very helpful in the search for the energy minimizers. 

Let us discuss the properties of functional (JHj). First, its extremals are stationary points 
of the corresponding nonlinear equation, equation 0, similar as in the case of the Rayleigh 
quotient (J2J) and equation (JBJ). (The unique stationary point is selected by the value of the 
/ 2 -norm, appearing explicitly in the energy functional as a parameter) . This follows from a 
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simple calculation: 
511 = 



d"*{ 6E %Psr + 6E Y f n 6f 



I 6f 



5^* 



ill* 



5f 



d n x' {ipSip* + ip*5ip) 



(10) 



where we have used that 



Sf 



d n x (if>6if>* + 



2||^|| 2 

By interchanging the order of integration in the double integral in formula (fTU|) and gathering 
the terms with 5tp* we obtain the variational derivative 



From equation (fTT|) it is quite clear that the stationary points of equation are extremals 
of the Rayleigh functional and vice versa. The corresponding chemical potential is equal to 
the integral in the second term on the r.h.s. of equation (jllj) . i.e. 



, , 8E{f, /*} 
V- = Re ( / d x — jji — / 



(12) 



which fact can be verified by direct integration of equation ®, i.e. fxf = 6E \^1 - , and taking 
into account the normalization of /. Equation (|T2|) plays here the role of the expression for 
the energy level e in the above problem of a quantum particle. 

The Rayleigh functional also distinguishes the local minima among the stationary points. 
To see this let us compute its second variation. Assuming that if) is a stationary point, we 
get 



sE{f,n 52f 



t— 

J -\\*\\ 



5 2 n^r} = (s 2 E{fj% 2f=0 + J d n x { 5E{ 6 f r n s 2 r + 

5 2 E{f,f% 2f=Q + J d n x {iif5 2 r + ^r5 2 f}^ 

www 

where we have used the fact that / satisfies the stationary equation [if = ' M ^ ■ Integrat- 
ing by parts and taking into account the normalization of /, i.e. 5||/|| 2 = 5 J d n x|/| 2 = 0, 



we arrive at 



5 2 K{i>,r} =[S 2 E{f,f% 2f=0 -2 f i I d n x\6f\ 2 



(13) 



IV-II 

The subscript in the first term on the r.h.s. of this equation means that the variable of 
variation of the energy functional is, in fact, f(x ). Taking this into account, one immediately 
recognizes on the r.h.s. of equation (j!3|) the second variation of the functional /*} 
defined in equation (JHJ) and evaluated in the space of normalized functions. Therefore, the 
minima of the Lagrange modified energy functional (for a fixed nonzero ^-norm) are also 
minima of the Rayleigh functional and vice versa. Importantly, the Rayleigh functional 
does not contain the trivial solution i/j — among its extremals, in contrast to the Lagrange 
modified energy functional. 

Finally, let us discuss the functional of Ref. which was proposed as another possible 
generalization of the Rayleigh quotient to the equations of the nonlinear Schrodinger type. 
In the latter work in the numerical search for the stationary solutions of two-dimensional 
Gross-Pitaevskii equation the following functional was employed 

FlM= I ^-*-J+? + m '» . d4) 
This functional can be expressed as follows 

^fl^WVH^"^ (15) 



Here the Rayleigh functional is given by 71 = J dx 2 f*(— V 2 +x 2 +NU\f\ 2 )f with / = VVIWI 
and N is the number of atoms corresponding to the stationary state. Note that the first 
variation of the Rayleigh functional is zero on a stationary solution and the last two terms 
in formula ([EH) have nonzero variation unless the norm ||^|| is kept fixed (H^H 2 = -^0- 
Hence, the functional ^{"0, ip*} introduced in Ref. p| cannot be employed for unconstrained 
minimization. 

Now let us show that, in fact, the Euler-Lagrange equation for the functional F{ip,ifj*} 
defined by equation (|14|) is self-contradictory (in other words, the first variation of F{ip, ip*} 
is never zero) and, hence, does not lead to any stationary solutions at all. The Euler- 
Lagrange equation for the functional (fTlj) reads 

(-V 2 + x 2 + 2U\ij\ 2 )ij = (A + F{^*})^, 



where F{ip, if)*} is the value of the functional on the solution. On the other hand, by 
multiplying the above equation and integrating we get 

(A + F{^*})|M| 2 = J dfV(-V 2 + f 2 + 2[/|^| 2 )^- 
hence 2U = U and we have arrived at a contradiction unless U — 0. Q.E.D. 

B. Generalization to several order parameters 

The Rayleigh functional was introduced above for a nonlinear system described by a 
single (complex- valued) order parameter if>(x). It can be easily generalized to the case several 
order parameters. The definition of the Rayleigh functional will depend on the number of the 
conserved Z 2 -norms. This number is determined by the type of phase invariance of the energy 
functional. We will concentrate mainly on the two broad cases: the incoherent and coherent 
coupling of the order parameters. For simplicity, consider the case of two order parameters: 
if> = (if)i,if>2)- In the case of incoherent coupling, there are two independent constraints 
corresponding to the two conserved ^-norms (the normalization is defined independently 
for each order parameter: ||/;|| = 1, I = 1,2). On the other hand, if the coupling is 
coherent, then there is only one constraint corresponding to conservation of the total Z 2 -norm 
(ll/i|| 2 + \\h\\ 2 = 1)- We also give an example of application of the Rayleigh functional to 
a nonlinear system which does not belong to either of the above cases (see section ITTTj) . 

An example of the incoherent coupling is the Gross-Pitaevskii functional for a two-species 
mixture of degenerate quantum gases in an external trap below the condensation tempera- 
ture (see, for instance, Ref. [l|): 

E^u^lMl} = I d 3 f |^iV i (|V^| 2 + A 2 x 2 |^| 2 )+ Yl ^f^NmW^m? 

U=l,2 l,m=l 

(16) 

On the other hand, the coupled-mode approximation for a Bose-Einstein condensate in the 
three dimensional parabolic trap modified in one dimension by a laser beam with creation 
of a central barrier illustrates the coherent coupling: 



E e {Mi,^,^} = n I d 2 x { V ( |W/| 2 + x 2 \H 2 + ^l^l 4 ) - KftMS + ^i) 



/d 2 xi W|V^| 2 + 

U=l,2 ^ 



(17) 
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for derivation of the coupled-mode system and the applicability conditions consult Ref. 

^]). The coupling coefficient k is proportional to the tunnelling rate through the central 
barrier of the resulting double-well trap. It is easy to see that functional (jl6|) admits two 
independent phase invariance transformations, ipi — ► e l ipi, while functional (|17|) admits 
only the simultaneous phase invariance transformation of both variables with the same 9. 
The Rayleigh functionals for the two cases corresponding to equations (|16|) and (jl7j) are 
defined as follows: 

v f jj^jtLJ^jtLX v ^ H iL A ill 

where ||^|| 2 = / d n x (|Vi| 2 + |^ 2 | 2 ). 

It should be stressed that to define the Rayleigh functional one should employ the nor- 
malization which follows from the phase invariance in its general form. For instance, the 
normalization by the total ^-norm in the incoherent coupling case would allow one to find 
only the stationary solutions with the same chemical potential for each component of the 
order parameter, i.e. a very limited class of solutions. 

The variational derivative of the Rayleigh functional TZ nc (fTSj) is derived by a mere repe- 
tition of the steps which led to equation (jllj) . We have 



8TZ nr . I f 5E nc f f 3 _>5E nc 



i mi 1 1 w 



(19) 



IIV'lll ' 11^211/ 

where, as usual, the /-variables are substituted in the energy functional instead of ?/>-ones, i.e. 
E nc = E nc {fi, /*, f 2 , /|}. Now it is easy to see that, thanks to the equation = 5E nc /Sf^, 
the integral on the r.h.s. of formula (|19p coincides with the corresponding chemical potential 
Hi. Thereby the extremals of the Rayleigh functional TZ nc are stationary points of the 
corresponding nonlinear system of equations. 

The variational derivative of the functional TZ C (|18j) can be most easily obtained from 
formula (jll)) by noticing that this Rayleigh functional depends on the function ip*(x ) through 
/* = and also through the total norm ||^|| in all other /-variables. Therefore, each 

variable of the energy functional leads to an integral term similar to the second term in 
equation (JTT|) but multiplied by ft. Hence the final result: 



5TZ C 1 SE, 



c 



<ty? \\n i we 




(20) 



1 11*1 
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where E c = E c {fi, /*, / 2 , / 2 *}- To show that the extremals of the Rayleigh functional 1Z C are 
also stationary points one has to use both stationary equations, i.e. /ifi = 5E nc /5f{ and 
pL j 2 = 5E nc /5f2, multiplied by the complex conjugate functions f*. We get 

/ d2f E 5 ^r m = »i^i\h\ 2 + \h\ 2 ) = ». 

J m=l,2 Jm ^ 

Let us now find the second variation of the Rayleigh functional at a stationary point in 
the above two cases. The derivation is similar to that for a single order parameter: one has 
to use the stationary equations and the normalization constraints for the /-variables. We 
get: 



5 2 U nc = (5 2 E nc \ s2f=Q -2 t*m I d 3 x|5/„ 

\ m=l,2 J 



?—( ^1 v>2 ~\ 

J -\\\i>i\\>M2\\) 



5 2 n c = ( s 2 E c \ s2f=0 - 2/i / d 2 x \ 6 U 2 

m=l,2 



(21) 



(22) 



F=-i- 

J \w\ 



In both cases the second variation of the Rayleigh functional coincides with that of the La- 
grange modified energy functional considered in the space of normalized functions. There- 
fore, the minima of the Lagrange modified energy functional (for the fixed values of the 
conserved Z 2 - n o rm s) are also minima of the Rayleigh functional and vice versa. 

Concluding this section we note that by introducing the concept of Rayleigh functional 
we have converted the task of finding the stationary points to an unconstrained optimization 
problem which also distinguishes between the local minima and the saddle points. In the 
next section we give examples of application of this optimization method to the numerical 
search for the stationary states of Bose-Einstein condesates. 



III. NUMERICAL SEARCH FOR STATIONARY POINTS 

By analogy with the use of the Rayleigh quotient in the numerical search for eigen- 
functions, the nonlinear optimization can be based on the Rayleigh functional. Leaving a 
mathematically rigorous analysis for future research, let us nevertheless give an heuristic ar- 
gument. The Rayleigh functional is a compound functional: the energy functional depends 
on f(x ) and f*(x), while these are functionals of ip(x ) and vp*(x). The numerical nonlinear 
optimization of a functional reduces to that of a smooth function (though of a large number 
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of variables). But a continuous function on a finite-dimensional sphere (which in the numer- 
ics plays the role of the variable f(x)) always has a minimum. Since the second variation 
of the Rayleigh functional is actually the second variation of the Lagrange modified energy 
functional, evaluated in the space of normalized functions, the optimization method based 
on the Rayleigh functional has the capacity to return the absolute minimum. Generally, it 
converges to a local minimum . 

An additional argument in favor of the nonlinear optimization based on the Rayleigh 
functional is provided already by the imaginary time relaxation method - a widely used 
approach in the computational physics. Indeed, it is precisely the Rayleigh functional which 
is minimized in this method. Let us recall the usual formulation given by equation (J1J of 
section |1] First of all, this formulation does not reflect the total change of ip, since after 
each numerical step one performs the direct normalization of the updated function. This 
means that the variational derivative 5E/Sf* does not tend to zero, since the norm of the 
order parameter is not forced to approach a constant value as r — > oo. However, one can 
reformulate the relaxation method for the variable f(x,r) = ip(x , r)/| \ip(x, r)| | which has 
a fixed /2-norm. We have 

df d f ip \ 1 fdip ip f ( dip* dip 

a x [ ip— h ip — 



dr dr\M\\) \\ip\\ {dr 2||^|| 2 ./ V dr T dr 



Sip* (23) 

I i> 1 1 

where we have used equation (£Q) to substitute for dip /dr. We see that if the imaginary time 
relaxation converges, i.e. if df /dr — > 0, then it converges necessarily to an extremal of the 
Rayleigh functional. 

Equation (|2*H|) . besides revealing a new side of the imaginary time evolution method, 
manifests its strong drawback: the order parameter evolves continuously and, hence, quite 
slowly. But, since it is actually a way to minimize the Rayleigh functional, a discrete 
minimization scheme can be adopted with the order parameter performing a series of finite 
steps to a local minimum. We have naturally arrived at the steepest descent formulation: 

^> = V>(« + D w = _L*nf ) ,f k) } (24) 

«fc 5ip*y k > 

where the parameter at is selected by some line-search algorithm. More advanced gradient 
methods, for example, the conjugate gradient method, can be employed for the finite-step 
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minimization. However, we have found that the steepest descent algorithm performs very 
well in the nonlinear optimization based on the Rayleigh functional. 

In selecting the step size l/a& we have adopted the Barzilai-Borwein two-point method 



3 step si; 



111 ] (see also Ref. [12|,ll3j]). There are two closely related Barzilai-Borwein methods: 

Re f d n x D*( k ~^ ( - f^j) 



f d n x 

(2) J 



Sip* W 



(25b) 



Ite/d»afZH*-i) ^^-^^ 

Both methods are fairly equivalent in terms of performance. The generalization of the 
numerical scheme fl24|) - (j2~5*|) to several order parameters is given by the substitution of the 
scalar function %p(x) by a vectorial one ip(x) defined similar to the vectorial function f(x) 
used in formulae (|T9"j) and (^Dj) (in the vector case, the inner product in equations (|25a|) and 



(I25bj) contains, besides the integral, the Hermitian inner product in the finite-dimensional 
target vector space). 

In the numerical implementation we have used the s pec tral collocation method for the 
spatial grid based on the Fourier modes (see Refs. j3QJ3)- The accuracy of the obtained 
solution was checked by direct numerical substitution into the governing equations. In all 
cases a rapid convergence to the stationary point was observed. We have cut off the allowed 
values of 1/at by a small positive number (in our case 1/a > 10~ 8 ) to have the iteration step 
always in the direction of decrease of the Rayleigh functional. 



Examples of numerical search for the stationary states 

A. Incoherent coupling. Our first example of the numerical optimization based on the 
Rayleigh functional is the calculation of the ground state in the two-species Bose-Einstein 
condensate of Rubidium in the quasi two-dimensional geometry (the pancake trap). The 
derivation of the governing equations and further discussion can be found in Ref. jlo| . 
The energy functional for the two-species Bose-Einstein condensate is given by the two- 
dimensional version of equation Ijlfijl . The parameters Ai and A2 account for the different 
Lande magnetic factors. For the two isotopes of Rubidium we have: = 8/7 (for 85 Rb) 
and Ag = 6/7 (for 87 Rb), while the interaction coefficients are gn = —0.0219, #22 = 0.0068, 
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and gu = 0.012. The variational derivative is given by equation f)19|) . To simplify the 
task, we note that the first term of the variational derivatives is the r.h.s of the stationary 
Euler-Lagrange equations written for the normalized functions: 

Hifi = N, (-V 2 /i + Xlx 2 h) + {gnNUM 2 + g 12 N 1 N 2 \f 2 \ 2 ) f u (26a) 
/i 2 / 2 = N 2 (-V 2 / 2 + A 2 f 2 f 2 ) + (feiVWaljil 2 + <? 22 A||/ 2 | 2 ) f 2 . (26b) 

Here the chemical potentials fij are defined as follows: fij — fjLj^.Nj, j = 1, 2, with fi N . being 
the chemical potential for the order parameter 0j normalized to the number of atoms Nj. 

A nontrivial feature of the two-species condensate is that the ground state suffers from 
the symmetry breaking transformation, if for instance, for a fixed number of atoms of the 
85 Rb isotope the number of atoms of the 87 Rb isotope increases In figure^we show the 
symmetry preserved (axially symmetric) ground states, while the symmetry-broken one is 
illustrated in figure 121 The axially symmetric states were found by the optimization method 
projected on the space of axially symmetric functions and formulated in polar coordinates. 
The polar radius grid contained 128 points, while the two-dimensional grid was 64 x 64 (we 
have used the 32 x 32-grid to produce figure EJ. The initial condition in both cases was the 
Gaussian profile (modulated by a symmetry-breaking perturbation in case of figure |2J). Our 
numerical simulations were performed in the MATLAB. The typical number of the iterations 
to reach the 10~ 10 convergence was between 400-1000 (on a personal computer with an AMD 
lGhz-processor the iterations have taken up to 5 seconds in the axially-symmetric case and 
up to 30 seconds in the 2D case). 

B. Coherent coupling. Consider the following functional 



^{^1,^,-02,-05} = / d 2 x {|V^| 2 + |W 2 | 2 + x 2 (|Vi| 2 + \H 2 ) +e\^ 



|2 



2 



-f l^ll 4 + - + 020!)}, (27) 

which appears in the description of the stationary states of a pair of repulsive and attractive 
(a > 0) two-dimensional condensates trapped in an asymmetric double-well potential ^] (in 
this case is proportional to the actual number of condensate atoms with the coefficient of 
proportionality of the order 10 2 — 10 3 ). Here e is the zero-point energy difference between 
the two wells of the external trap and N is the total number of atoms. The Euler-Lagrange 
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equations for the energy functional (}2Tj) read: 

/i/i = -V 2 /i + x 2 h - Kf 2 - iV|A| 2 /i, (28a) 
Hh = -V 2 / 2 + ef 2 + x 2 f 2 - Kh + aN\f 2 \ 2 f 2 , (28b) 

where /i coincides with the chemical potential for the order parameter normalized to the 
number of atoms N. 

As in the previous example of the incoherent coupling, we have observed a rapid con- 
vergence of the iterations to the ground state solution. The numerical minimization was 
performed in polar coordinates. The characteristic times to reach the deviation of about 
10~ 10 are the same as in the previous example (formulated in polar coordinates). To verify 
that the solutions obtained with the use of the Rayleigh functional method are indeed the 
ground states, we have used another numerical approach which enables one to find all possi- 
ble stationary states (consult for more details Refs. jslH]). The deformation of the ground 
state with the variation of the zero-point energy difference is illustrated in figure El where 
the condensate of vanishing atomic interaction, a = 0, is tunnel-coupled to an attractive 
condensate. 

C. Arbitrary coupling. The minimization method based on the Rayleigh functional can be 
applied to systems of equations coupled by a different way than the coherent and incoherent 
coupling. As an example, we consider here the Bose-Einstein condensate in a pancake trap 
(in the r± = (x, y) plane), but in contrast to the above considered cases, we take into account 
the transverse dimension (z) by using the Gaussian variational ansatz with the r*j_-dependent 
parameters, 

V> = T74 iL~ a ex P f- n S /fa. ty aM2t \ (29) 



A similar ansatz was already 



^/4 w i/2(f 1)t ) ^ 2w 2 (r ± ,t), 
for the order parameter in the Gross-Pitaevskii functional 

used with success for the condensate in a cigar shaped trap |l9j. Here the parameter az 2 /2 
accounts for the transverse motion of the condensate. We are interested in the stationary 
solutions which correspond to a = and are given as / = e-^^U^) and w = w(r_j_). 
Integrating over the transverse coordinate in the Gross-Pitaevskii functional (J2J) we get the 
energy functional for the stationary states: 



E 2D = j d 2 r ± |^-|V^| 2 + y(rl)|C/p 



,2 



2v / 27r w Am \ w 2 w 2 a\ 

(30) 
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Here a z = {h/muj z ) 1 / 2 the transverse oscillator length, V(r±) = mu 2 _r]_/2 is the external 
potential in 2D. Note that w — > a z as |r_j_| — > oo. The pancake geometry corresponds to 
A = uj z /ujx_ ^> 1. The equations satisfied by the stationary states read: 

,u + fv{u - v(f L) u - >\ir?u - £ (S^ + \ + ^' 

2m V27TW 4m V ur ur a* 



^J^\uA + + -1 - ^) |I7|» + = 0, (32) 

where a s is the s-wave scattering length. 

Assuming that the transverse width w is of the order of a z (which is supported by the 
numerics) and estimating the operator Vj_ ~ l/a± we see that the term U(V ±w) 2 /(w 2 ) is 
much less than (^ + ^)C7 for A 3> 1. If this small term is omitted from the energy functional 
((HOI), the equations for the stationary state become equivalent to the two-dimensional non- 
polynomial NLS equation (2D NPSE) of Ref. H- 

(To compare with the latter reference, 
though, one should take into account that the coefficient 1/2 is missing there at (^ + 
in the equation which is equivalent to equation (|3Tj)). 

The essential difference of the energy functional given by equation (}30|) from the above 
considered energy functionals lies in the fact that it has two variables of variation, U and w, 
but only one phase invariance exists and is associated with the normalization of U. However, 
the minimization method based on the Rayleigh functional defined with the substitution 
/ — ► £//||[/|| into the energy functional (|3Uj) works in this case also. To explain this, we 
note that one can think of the "chemical potential" corresponding to the variable w as being 
equal to zero, thus the variation of the Rayleigh functional with respect to w coincides with 
that of the energy functional. 

The solutions of equations (pTTj) and (jS2j) were searched for in polar coordinates. As the 
system similar to (|31 )) -([32 )) was proposed in Ref. j3] as an improvement to the NLS equation 
in two dimensions, we aimed to find out when does one need to use this complicated system 
for the pancake condensate instead of the much simpler 2D NLS equation? We have found 
that only in the study of collapse instability for g < the system (|31 jl -(|H2 jl gives visible 
improvement over the 2D NLS equation. In all other cases, if one is not interested in the 
transverse profile of the condensate, the 2D NLS equation gives satisfactory results for the 
stationary states. (This, however, may not be so for the arbitrary time-dependent evolution, 
which is not captured by the 2D NPSE system of Ref. (2(| since their a = identically.) 
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To illustrate this, in fig. Hjwe give the relative difference between the solutions U(r±) 
of the system (j31j) - (}3*2*j) and the corresponding 2D NLS equation for several values of the 
nonlinear ity strength in the case of repulsive condensate (g > 0). The solution is given in 
the units of l/a±, a± = (h/muj^) 1 ^ 2 , the radius r± in terms of a±, and fi in terms of Jtlu±/2. 
The dimensionless nonlinearity strength is defined as G = 4^/2na s /a z . We have set A = 100 
which corresponds to the ratio a±/a z = 10. The transverse width w, given in terms of a z , is 
shown in fig. The maximal difference from a z is of the order of 10%. 

Though in the case of attractive condensate the relative difference between the solutions 
of the system (pTTf - ()32j) and the 2D NLS equation is also small, the critical value for collapse 
differs noticeably. The critical value of the number of atoms is defined by dfi/dN = oo, i. e. 
it is the endpoint coordinate of the curves in fig. H3 The system (|3~T]) - (j32j) is an improvement 
over the NLS equation in this case, since the critical value differs by only 2% for A > 20 



2l|. The 2D NLS value is 



from the actual value of the 3D condensate obtained in Ref. 
achieved as the asymptotic limit when A — > oo. 

We have argued that the minimization based on the Rayleigh functional is always suitable 
for obtaining the ground state solutions. However, the possibility of obtaining excited states 
was not excluded at all. In contrast to the linear eigenvalue problem of section |H] in the 
nonlinear case there is no general approach which would guarantee computation of all excited 
states, though there are some advances in this direction [5j. However, a proper selection of 
the initial condition for minimization may lead to convergence to the "nearest" excited state. 
We illustrate this on the multi-vortex solutions of a nonlocal Gross-Pitaevskii equation, 



which in dimensionless 
(see the details in Ref. 



orm, in the frame of reference rotating with the frequency Q, reads 



22j) 



id t m = -V 2 ^ + x 2 ^ - ttL z ^ + 9^ J dx' 2 K(\x - x'\)\V(x ')\ 2 . (33) 

Here K = -^iKq (^), Kq(z) is the Macdonald function, a is the effective range of atomic 
interaction, L z is the angular momentum projection on the z-axis: L z = —i(xd y — yd x ) = 
—idg, where 9 is the polar angle, and Q is the rotation frequency. For simplicity, we consider 
the periodic boundary conditions along the z-axis, ^(z + d) = ^(z). 

We searched for new non-axial vortex solutions, involving combinations of vortices and 
antivortices, which represent the excited states of the system, since in the nonlocal GP equa- 
tion ()33|) . similar as the in local one, the stationary vortex solutions involving combinations 
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of vortices and antivortices are never the energy minimizers. The energy is minimized by 



221. 



Tkachenko lattices comprised of vortices only 

In figures [7| and |S1 we give two examples of excited states comprised of vortex and anti- 
vortex combinations (the right panel) and the respective energy minimizers (the left panel). 
These solutions were found by using the initial conditions with the phase resembling that of 
the multi- vortex solution in quest. 



IV. RAYLEIGH FUNCTIONAL IN THE ANALYTICAL APPROACH 

The ground state solution of a nonlinear partial differential equation is degenerate if fixing 
the /2-norm does not select a unique chemical potential. In the critical case, when there are 
infinitely many stationary points with equal energy and the same Z 2 - n orm, the nonlinear 
optimization method based on the Rayleigh functional may not converge at all. 

To clarify this let us consider the critical NLS equation. The ground state of the focusing 
critical NLS equation is infinitely degenerate. It is well-known (see, for instance, Refs. 

that the two-dimensional critical NLS equation has a family of soliton solutions 
having the same ^-norm HV'II 2 ~ H-69 and a specific member if) = e lt R(x) satisfying the 
stationary equation 

V 2 R + R 3 -R = (34) 
is called the Townes soliton. The stationary NLS equation 

fiif) + V 2 ^ + \il)\ 2 i) = (35) 

admits the scale invariance given by ip(x) — > kip(kx) and /i — > k 2 fi. Easy to see that, 
in two spatial dimensions, the scale invariance preserves the /2-norm of the solution, while 
the energy of the stationary solution is exactly zero (see below). This is precisely why 
one cannot obtain the Townes soliton by the nonlinear optimization method based on the 
Rayleigh functional - the method fails to converge. 

This drawback of the Rayleigh functional in the numerical approach turns out into an 
advantage for the analytical analysis. Indeed, for the critical NLS equation it immediately 
leads to the well-known identity [18] 

~H J d 2 x |V>| 2 = J d 2 x |VVf = \ J d 2 x\tfj\\ (36) 
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which is satisfied by any stationary solution. Let us show this. First of all, it is easy 
to see that the energy functional of the NLS equation E = J d 2 a;(|V^| 2 — |^| 4 /2) has a 
scale invariance property: E — > k 2 E for ip(t,x) —* kip(k 2 t, kx). The invariance property is 
transferred to the Rayleigh functional: 1Z — > k 2 lZ. Differentiation of the latter with respect 
to the scale invariance factor k gives 1Z = at the stationary point. Hence, the second 
equality in equation (jHUj) follows. The first equality is derived by integrating the stationary 
equation multiplied by tp* and using the already proven equality. 

In the above derivation we have explicitly related identity (J36J) to the scale invariance 
property. The scale invariance is responsible for the exact balance of energies in any sta- 
tionary solution of the critical NLS equation: the kinetic energy is balanced by the energy 
due to the self-interaction. 

The identity (j3T?j) is well-known. However, the point we try to make here that its relation 
to the critical scaling makes such an identity universal, i.e. it appears in connection to any 
nonlinear equation for which the energy functional has a scaling invariance for a family of 
solutions with the same /2-norm. Such is also the one-dimensional critical NLS equation (in 
the stationary form) 



E — > k 2 E, with E = J dx(\difj/dx\ 2 — |?/>| 6 /3), while the ^-norm N — J dx\ip\ 2 remains 
unchanged. Hence the Rayleigh functional scales as follows 7Z — > k 2 lZ and, by repetition of 
the above arguments, we get a similar identity satisfied by the family of stationary solutions: 



Note that one cannot use the energy functional in the above argument: the stationary 
solution is not an extremal of the energy functional. Using the Lagrange modified functional 
will not help either due to the explicit appearance of fi. The use of the Rayleigh functional 
is indispensable in this short derivation. 





(38) 



V. CONCLUSION 



Numerical search for stationary points of nonlinear partial differential equations is a 
difficult problem. To tackle such a problem one is left to try various methods and choose 
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the one which gives a better performance and accuracy. In this paper we have proposed a 
new numerical method, the nonlinear optimization method based on the Rayleigh functional. 
This method is a natural generalization of Poincare's minmax principle for linear equations, 
formulated with the use of the Rayleigh quotient. 

It turns out that the imaginary time relaxation method, a widely used method in the 
computational physics, is just a special case of nonlinear optimization based on the Rayleigh 
functional. We have used the gradient scheme for minimization of the Rayleigh functional, 
but this is not essential: one can use the minimization schemes which do not require use of 
the gradient. 

The simplicity of the Rayleigh functional and its universality for nonlinear equations is 
one of the advantages of the method. Moreover, the method can distinguish between the 
local minima and saddle points, since the second variation of the Rayleigh functional is 
equal to the second variation of the Lagrange modified energy functional, if the latter is 
evaluated in the space of normalized functions. The Rayleigh functional takes care of the 
normalization constraint of the stationary solution. Other constraints, however, have to be 
treated in the usual way. 

There are, however, some exceptional cases when the nonlinear optimization based on 
the Rayleigh functional may fail. The principal cause of the failure is the infinite degeneracy 
of the ground state solution. Still, a failure in the numerics is partially compensated by the 
fact that the very degeneracy allows one to get an important analytical insight relating the 
scale invariance and an identity satisfied by the stationary solutions in the critical case. 
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VI. FIGURES 
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FIG. 1: The axially symmetric ground state of the two- isotope Bose-Einstein condensate mixture 
of Rubidium. The dashed and dotted lines give the order parameters (normalized to the actual 
number of atoms) of 85 Rb and 87 Rb isotopes, respectively. Panel (a) = Ngf = 300, /U85 = 1.18 
and pgr = 2.67; panel (b) N 85 = 200, N S7 = 15000, /x 85 = 11.69 and /i 87 = 7.827; panel (c) 
iVgs = iVg7 = 500, [i85 = —3.48 and /X87 = 3.08 (in dimensionless units). 
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FIG. 2: The symmetry breaking ground state of the Bose-Einstein condensate mixture of the 
isotopes of Rubidium. The left and right panels give the shapes of the order parameters of 85 Rb 
and 87 Rb isotopes, respectively. Here = 150, Ngj = 20000, fi^5 = 13.19 and /j,s7 = 8.91 (in 
dimensionless units). 




FIG. 3: Deformation of the ground state of two Bose-Einstein condensates trapped in an asymmet- 
ric double-well potential with variation of the zero-point energy difference. The picture corresponds 
to the energy functional (|27|). Here a = 0, k = 1 and N = 15, the solid line gives ipi, while the 
dashed one tp2- 
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FIG. 4: The relative difference between the solutions of the 2D NLS equation and the system 
l[3T ]l -(|33 |l . Here A = 100. 
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FIG. 5: The transverse width of the condensate given by the system (|!TT j) - p2[) , Here A = 100. 
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FIG. 6: The dependence of the chemical potential on nonlinearity for the ground state of the 
system (|31 j) -([32 j) . The dashed line gives the same for the 2D NLS equation with the parabolic 
potential. 
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FIG. 7: The stationary vortex solutions for a = 0.052, g = 10 4 and f2 = 0.38. The left panel 
shows the 3- vortex solution and the right one - the 5- vortex solution comprised of 4 vortices and 
1 antivortex (in the center). 
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FIG. 8: The stationary vortex solutions for a = 0.052, g = 10 and O = 0.5. The left panel shows 
the 8-vortex solution and the right one - the 16-vortex solution comprised of 12 vortices and 4 
antivortices. 
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